{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "40e71b28",
   "metadata": {},
   "source": [
    "# 期末项目 Final Project\n",
    "\n",
    "- 姓名：杜东洋\n",
    "- 学号：20201050438\n",
    "- 邮箱：Starry_339@163.com"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "418b1850",
   "metadata": {},
   "source": [
    "$$\n",
    "n(x) = (\\frac{1}{2}(2x)^{\\frac{1}{2}} * sigma)  e^{\\frac{- x ^ 2}{2sigma ^ 2}}\n",
    "$$\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "7a19880b",
   "metadata": {},
   "source": [
    "## 一、导入数据库"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "a9b448bc",
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "3401f7b6",
   "metadata": {},
   "source": [
    "## 二、代入参数"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "1f819622",
   "metadata": {},
   "outputs": [],
   "source": [
    "D = 1\n",
    "M = 40 * D\n",
    "E = 10\n",
    "tau = 0.01\n",
    "dx = 0.001\n",
    "dt = 0.00001"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "45b6caca",
   "metadata": {},
   "outputs": [],
   "source": [
    "sigma = D * dt / dx ** 2\n",
    "gamma = M * E * dt / dx"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "1fed12ad",
   "metadata": {},
   "outputs": [],
   "source": [
    "a = - sigma / 2 \n",
    "b = 1 + sigma + dt / tau + gamma\n",
    "c = - sigma / 2 - gamma\n",
    "d = -a\n",
    "e = 1 - sigma\n",
    "f = d"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "2c3505e5",
   "metadata": {},
   "source": [
    "## 三、构建方程并作图"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "fc1ac54a",
   "metadata": {},
   "outputs": [],
   "source": [
    "m = 500\n",
    "xBegin = - 2 * m * dx\n",
    "xEnd = (2 * m + 1)* dx\n",
    "xArr = np.arange(xBegin, xEnd, dx)\n",
    "\n",
    "n = 5000\n",
    "tBegin = 0\n",
    "tEnd = dt * n + dt\n",
    "tArr = np.arange(tBegin, tEnd, dt)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "b88bc122",
   "metadata": {},
   "outputs": [],
   "source": [
    "width = 10 * dx\n",
    "uBegin = 1 / np.sqrt(2 * np.pi) / width * np.exp(- xArr ** 2 / 2 / width ** 2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "9840646a",
   "metadata": {},
   "outputs": [],
   "source": [
    "A = np.identity(xArr.size - 2)\n",
    "B = np.roll(A, 1)\n",
    "C = np.roll(A, -1)\n",
    "B[0, 0] = 0\n",
    "C[-1, -1] = 0"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "id": "8ab1142f",
   "metadata": {},
   "outputs": [],
   "source": [
    "A1 = b * A + c * B + a * C\n",
    "A2 = e * A + f * B + d * C\n",
    "A3 = np.linalg.inv(A1) @ A2"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "id": "4b2c75c0",
   "metadata": {},
   "outputs": [],
   "source": [
    "uArr = np.zeros((xArr.size - 2, tArr.size))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "id": "0551e606",
   "metadata": {},
   "outputs": [],
   "source": [
    "uArr[:, 0] = uBegin[1:-1]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "id": "ddebe115",
   "metadata": {},
   "outputs": [],
   "source": [
    "uCurrent = uArr[:, 0]\n",
    "for i in range(1, n):\n",
    "    uNext = A3 @ uCurrent\n",
    "    uArr[:, i] = uNext\n",
    "    uCurrent = uNext"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "id": "6a57d9f1",
   "metadata": {},
   "outputs": [],
   "source": [
    "XX, YY = np.meshgrid(tArr, xArr[1:-1])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "id": "1719f7b7",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(2000, 5001)"
      ]
     },
     "execution_count": 13,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "XX.shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "id": "4ecfabdb",
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(2000, 5001)"
      ]
     },
     "execution_count": 14,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "uArr.shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "id": "41ffc08e",
   "metadata": {},
   "outputs": [],
   "source": [
    "uArr = np.zeros((5000, 5000))\n",
    "for i in range(5000):\n",
    "    for j in range(5000):\n",
    "        uArr[i, j] = 255 - i + j"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "id": "81f9c975",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.colorbar.Colorbar at 0x161cd297490>"
      ]
     },
     "execution_count": 16,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAUwAAAD8CAYAAAAc052eAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/MnkTPAAAACXBIWXMAAAsTAAALEwEAmpwYAAAj7ElEQVR4nO2df8xk1XnfP8/78sOOHRJjDNnuuyo02VTFNLXDaruVq6o1pmxsi+WPWFpLLivV0koWkWy5lQtFqpQ/VrJbKbZQiitkR1lkNwQlqVhZAZdAUFSJgtc/ML+8ZjFp2EJZYfKDODJm33n6xz135tw7d2bOzJyZ9x7v97s6mnOf89znOfPrs+fc+8455u5IkiRJs7Wx0x2QJEkqRQKmJElSogRMSZKkRAmYkiRJiRIwJUmSEiVgSpIkJWrtwDSzg2Z2ysxOm9mt684vSZK0qGydf4dpZpvA94HrgTPAN4CPuvsza+uEJEnSglr3CHM/cNrdf+DuPwHuAQ6tuQ+SJEkL6YI159sNvBgdnwH+advJzI4CRwHsoouuvfCKy6dHtfhgyojZJjcl+diM+LXbrDzm6d3pyGkTDmxC36wjWZdv7DexnVF/zKa0dfXDRjHj9nYuiyrjbTNiWLvNx/rSHcMbeSf5NW2T3pe0WNNjTMrZVPWe2eS2Yd067QD/58VzvPradsrHcaJu+Fdv8x++tp3k+83vvvF1dz+4TL6d0rqB2fWmjH0z3f0u4C6Ai/fs8b/36U9VZ4azffhJ8vhTNXx08zHbqAcdbVNsZuPnWcvfYv/a3XwshkXHo/oUW4i1YT7WtmFxfdS20XFuO0Zs26Ad1xvtwNBnwwbjbVGMkW3A5ti509pGcTcZjMVv2zZtMDx30wbDWMN6eNyM4sf1Ok87ZxW3rrf7M4oxyjkY1ofxu2xRznGbRzZvxfKoH7UNNsObW08PN4HN8CEb2cKx2ageWjfNhvUNjAMHz7CsXn1tm8e+vpXke+Gu5y9bOuEOad1T8jPAnuh4C3gp6UxniNbhQMVthNvo0Wqittvqc9ptU2zu4+d5y99j/9rdbSyGR8ej+hRbiDVwG2sbeFwftQ06zm3HiG0D2nGt0Q4MfQa+Md4WxRjZNtgeO3da2yjudvhIxvHbtm3fGJ677RvDWMN6eNyO4sf1Ok87ZxW3rrf7M4oxyrkxrA/jd9minOM2i2zWimVRP2obbIc3t0IubAPb4UM2soVj91E9tG67D+uDhBlTmpxtHySVkrVuYH4D2GtmV5nZRcBh4MSsk6wBsrZN0Kzb6zZBU9CsbGnQzCGngm9KKVlrBaa7nwN+A/g68Cxwr7s/nXKuoCloCpr9huYg8V/JWvvfYbr7H7v7L7v7L7r7sbSTqgdBU9AUNPNDM8eYz3He9EFSKVnl/NJH0BQ0ETRXAc0c8pAjpZSsMoDZApKgKWgKmv2Dpq5h9kTWATBBU9AUNPNAM8ek3EOulFKy+g/MGIaCpqDZii9o5oFmDg0SS8nqPzBB0JxkC7EETUFzGWjmkCdev9Q1zHVJ0BQ0Bc2VQDMHwtzhzcRSsnoPTIMxqAmagqagmReay8uGfZxVSlbvgQmCpqApaNY+q4BmDjnV5yyllKwigAmCpqApaNY+fYWmRph9UQ2+qC5oCpqCZh5o5hj0OQJmvyRoCpqC5sqguawceNM3kkrJKqP3bUB22QTNpi3EEjQFzVnQJAM0HQt9nF1KVjm9FzQFzYltgmbj+c0JzTxjzNF7P6uUrP4D01vQix4FTUFT0MwDzWXlUT90DXOnJWg2YgiagmZOaOZR9bqllJJVTu8FzUYMQVPQ7BM0PfQ3pZSsMnrfANi4DQRNQVPQrH3mhWYOZrobP/HNpFKyygAmCJoImoLmaqCZSwMsqZSsIoDZBUhBU9Cc3SZoNp7fCqHpw37k+bMiM9s0s2+b2dfC8aVm9qCZPRce3xH53mZmp83slJndENmvNbMnQ9sdZrb0ky0CmLgJmnEIQbN1rqC5DDQ9CzSz3/T5JNVGibVuBR5y973AQ+EYM7uaavfZdwMHgTvNrJ73fxE4CuwN5eCyz7IMYIKgKWgKmpNsGaC5rDz0LcdNHzPbAj4EfCkyHwKOh/px4KbIfo+7v+HuLwCngf1mtgu4xN0fdXcH7o7OWVhFADMGmKApaNZ+gmYeaMYf3WW07ZZUEvQF4DPQ+DPRK9z9ZYDweHmw7wZejPzOBNvuUG/bl1L/gdkBMEFT0Kz9BM3loZlDjvGmX5BUgMvM7GRUjtZxzOzDwFl3/2Zi6q4n4VPsS+mCZQOsReHpmwc2uGF4qFdteHho2bBQbdu8epNDpRU/5G3lHCWJcxoexWi8Vd6MX8WdbHMHs+Z5jmGRv4ecmI9CuFUpoxju1Xke4lf1KbYQa+DGBjTaBg4b1HVjI7QNoKrH59Y+IUZsw2CDOG71Qm3EfhgbePiiD5ptUYyRrfLbbJw7rW0Ud5sNNhm04jdt274BNmCDqr5pgwqidd032LAB2xibdXwbDOtV2wZ4M2cVFzYYsO3GZqM/oxijnBvgsBnn6rLFOWn3o3q9N6PXKPdNn0S96u77JrS9D7jRzD4IvAW4xMy+ArxiZrvc/eUw3T4b/M8Ae6Lzt4CXgn2rw76U+j/ChImjPo00NdKs/TTSXHykmWNK7qRNx2dNyd39NnffcvcrqW7mPOzuHwNOAEeC2xHgvlA/ARw2s4vN7CqqmzuPh2n762Z2INwdvzk6Z2GVAUwQNAVNQZPVQTOHVvxLn88C15vZc8D14Rh3fxq4F3gGeAC4xd23wzmfoLpxdBp4Hrh/8WdXqYgp+aypsqbnmp5rer749DyH3EdQzyV3fwR4JNR/CFw3we8YcKzDfhK4Jmefihlhzhr1aaSpkWbtp5Hm/CPNZVXd9NlMKiWrDGAmAkzQFDRrP0EzHZp5/nA97y99+qpyei9oCpoImquAZg45aYsH57xmuhPqPTAtgpSg2fQRNAXNKmc/oKkRZk8kaAqagubqoBl/JBeVU792s0vJKqb3gqagKWiuDprLK217itxLyq1bZQCzAbAOW30saAqaCJrzQjOHHHSXvFcSNAVNQXMl0Mxxl7x6jzUlx8x+x8zOmtlTkW1nFvMUNAVNQXMl0MyhzOth9lIpvf9dxhfeXOtintaASmQTNAVNBM12jHmh2fj8LSgPzzGllKyZwHT3PwNea5kPsebFPAVNQVPQXCE0l5a22Z2mlS7maWZH67Xyzv3dj8YhBIJmyyZoCprtGKnQzDDArEaYrj9cn1ddr4ZPsXfK3e9y933uvu+Cn3lbw1vQFDQFzfzQXFb6Lfl0vRKm2ax8Mc8OkAmagqagmQ+a3eOZ+bXi5d16oUV7f4J1LuYpaI63CZqCJnmgGX/UFlW1vFu2PX16q5Q/K/o94FHgH5rZGTP7ODuxmKegOd4maAqa5IFmDp0P1zBnLiDs7h+d0LS+xTydilQeHmxkMx8xoWFjdE7Dr/5GNmxahBjTIsTn6yLEOVStVlT2dDtF5TzDKaM/jTQ10tRIc7mR5rJy4E3fSColq4jep4BM0BQ0Bc3FoJlnAWH9NLJXEjSjx/qcdpugKWgyPzRzSb/06Yu6oBMeBc1Wm6ApaLJ+aOoueY9knYAZPQqarTZBU9AkHZq59vTRlLwP6gCeoBk91ue02wRNQZN0aC4r7enTJwmagqaguTJoxh+jReXAOd9IKiWrnN4LmoKmoLkSaOaSpuQ90GT4RPWWj6DZahM0BU1WDE3XlLw3mgg3QVPQRNAc5lwQmlm2qAj9058V9USCZldOQVPQzAPNHNIIsy+aBTdBU9BE0BzmnBOajc/MgnIEzH5J0BQ0Bc3VQXNJOca5wUZSKVn97/0k4AmagiaCZg5oZhhgDvuja5h9kKA5llPQFDTrWMtCkxwQc03J+yVBcyynoClo1rGWh+ZycgTM/mgW8ARNQRNBc1Foxh+PZSRg9kmCpqApaK4MmsvKMbYHG0mlZBXRext+2xg+CpqCZmUTNOt6HWteaOZSrps+ZrbHzP7UzJ41s6fN7JPBfqmZPWhmz4XHd0Tn3GZmp83slJndENmvNbMnQ9sdYRPGhVUEMEHQTMkpaAqadax5oJnllz6edUp+Dvh37v6PgAPALWZ2NXAr8JC77wUeCseEtsPAu4GDwJ1mVm+A/kXgKNUOtntD+8IqA5hDOAiagqaguQpo5pCH13RWmR3HX3b3b4X668CzwG7gEHA8uB0Hbgr1Q8A97v6Gu79AtTPtfjPbBVzi7o+6uwN3R+cspCKA2QSeoCloCpo5oZlHaaPL0N/LzOxkVI5OjGp2JfBe4DHgCnd/GSqoApcHt93Ai9FpZ4Jtd6i37Qtr5ja7O64ICA5gFTTdAhis8jHCd9Frn5F/l8286T8xhkX9mBa3Ph6Lry18MW3h2+ctfHMpZfQY9Kq775vlZGZvB/4Q+JS7/82Uy49dDT7FvrCKGGF2jxI10tRIUyPNHCPNpQgS5A7bA0sqKTKzC6lg+VV3/6NgfiVMswmPZ4P9DLAnOn0LeCnYtzrsC6sMYIKgiaApaK4Gmt0DsfmV8S65AV8GnnX334qaTgBHQv0IcF9kP2xmF5vZVVQ3dx4P0/bXzexAiHlzdM5CKgKY0yAoaAqaIGguD83l5JDtpg/wPuDfAO83s++E8kHgs8D1ZvYccH04xt2fBu4FngEeAG5x9+0Q6xPAl6huBD0P3L/M8+z/NUyovuAw8Rqlrmnqmqaha5qLXtOM/69cXPl+xePu/4vJw97rJpxzDDjWYT8JXJOlYxQywgRmjhw10tRIEzTSXHSkmUPV+zK7lKwygJkIQUFT0ARBc15oZhpi5pyS91a9B+bw5RU0BU1Bc3XQXFLVXXL9lrwXSgaYoCloRrEEzTRo5vhpZPVapJWS1X9gzgswQVPQjGIJmrOhmUuakvdFgqagCYLmxLadh6aTBksBcw1aCGCCpqAZxRI0J0MzfmuXkSeWkjUTmH1Zm07QHLeBoClo5oHm0nLwgSWVkpUywtz5temWAZigKWhGsQTNcWjmkqbkQG/WphM0BU1BcyXQ1F3ydM3138u61qYzs6P1WnnnfvyjGRBstgmaXTkFTUFzOjSXlYfne96PMGu116ab5tph8yn2caP7Xe6+z933XfCWtyVAsBVN0OzIKWgKmt3QzKL6vUopBSvpFdvRtemSITivv6ApaAqaOaGpKTn0Ym26GG6C5vRcgqagOS8081zDTLtDfj7cJe/F2nSCJoKmoLkSaDbeu2XkiaVgzVwPsxdr03nVg3jdSyN8V4Zt4b2o61HbbH+tpzkr5yhJnFPraTbOtTLX08yi+D+5n2IV+EufeiihkeasXBppopFmwkgzfsuWkieWglUEMEHQFDRH/oJmfmjmkSWWclUGMGfATdCcnkvQFDSnQTObBomlYPUemMO3VNAUNAXN/kKzfi9SSsHqPTAhHW6C5vRcgqag2QVN/TQyXf0H5gT4CJptm6DZCCFots6d3JZNnlgKVv+BCYJml5+g2fAXNHsAzVlT8ej9KFVFAHMafATNtk3QbIQQNFvnjrfFb8syMk8rJasIYIKgKWgKmrA6aC4tNxgkloJVBjDDp1zQnOAnaDb8Bc35oJlNnlgKVhnABEFzWtz6WNAUNJkfmsPXblkJmP3Q6MsoaAqaguYqoJlFAmZP5IJmUtz6WNAUNEmHZhbVr3NKKVj9B2bnF1zQFDQFzb5B0zytlKz+A5Pml1LQTIhbHwuagiazoZlNnlgKVhHABEGz4ZMStz4WNAVN1gNNjTB7Ie/8UgqaCXHrY0FT0GQyNLMxrH4dZ5WCVQAwA4wETUFzGFfQzAlNclzH9DlKwSoCmCBoCpqtNkEzMzQzKBMwzeygmZ0ys9Nmdmu+Di6vMoAZw0jQFDSHcQXNHNBMYFiSbJBWpsYw2wT+K/BrwNXAR83s6kxdXFq9B+Y4JATNMX86bIKmoEk6NLPIE8t07QdOu/sP3P0nwD3AoXydXE69ByZE/ysJmoKmoJkdmjlknl5maDfwYnR8Jth6oSKACYJmUgw6bIKmoMl0aHr82i6j+vWaVeAyMzsZlaNRlC6C5+rh0pq5L3kv5IBV0PSN6Nidan/t6gvoVrUZdd3BKgg4jPkNj4niad9ztO/5+bXveTalY+1Vd983oe0MsCc63gJeWqJXWVXECLMxmtRIc3YMOmyT/Drja6TZCPFTPtLsHtTNr0xT8m8Ae83sKjO7CDgMnMjSwQwqApggaAqa4/GH57TbBM0FoLmknCx3yd39HPAbwNeBZ4F73f3pPJ1cXmUAswuMgqagSXROu03QTIZm/FIuJU8ss8K4/7G7/7K7/6K7H8vVvRzqPzDHvoCC5lw56bAJmoIm5P2TorpzGYDZZ/UfmCBoLpuTDpugKWiSF5qZrmH2WkUAs/sLLmjOlZMOm6ApaEpzqQhggqApaLZzCpq9g6YnloJVBjA74CZoCpqCZh5oeqbVinLcJe+7ygAmCJoTn/sCOemwCZrnNTSzyBNLweo/MGfATdAUNAXN5aCZA2JG9T6llJLVf2ACDFzQnOAvaNY5Bc1FoZmNYZ5YCtZMYJrZW8zscTN7wsyeNrPfDPZLzexBM3suPL4jOue2sPjnKTO7IbJfa2ZPhrY7zCx9PiBoTvQXNOucguai0FxarhFmrTeA97v7PwHeAxw0swPArcBD7r4XeCgcExb7PAy8GzgI3BkWBQX4InAU2BvKwZRO2uhTIGhO8Bc065yC5rzQzKZBYilYM4Hplf42HF4YilMt6nk82I8DN4X6IeAed3/D3V8ATgP7zWwXcIm7P+ruDtwdnTNTguZsf0Gzzilo7gQ0NcIMMrNNM/sOcBZ40N0fA65w95cBwuPlwX3SAqC7Q71t78p3tF4r7803fzT+RRE0J/oLmnVOQTMVmtnkiaVgJQHT3bfd/T1Ua9PtN7Nrprh3vQM+xd6V7y533+fu+y668G0NT0Fztr+gWecUNFOgmevvMAXMltz9r4BHqK49vhKm2YTHs8Ft0gKgZ0K9bZ+pNjAEzdn+gmadU9BMgWYOaUoOmNm7zOznQ/2twAeA71Et6nkkuB0B7gv1E8BhM7vYzK6iurnzeJi2v25mB8Ld8Zujc6bLBc123BR/QbPOKWhOg2Y2eWIpWClbVOwCjoc73RtUC3p+zcweBe41s48DfwF8BMDdnzaze4FngHPALe6+HWJ9Avhd4K3A/aFMVwwfmls9mIfP4MCxDavq7rS3pRj5R1taaLuL2THq79OsuPXxWHxtd4H1f7uLXCr9Z48pmglMd/8u8N4O+w+B6yaccwwYW/jT3U8C065/dqoJH0FT0BQ0c0Izy6jPyROn5yrjlz7QmuaO6lVbcNL0fKL/UjnpsE3y64yv6XkjRM+m5zk4Z3OUklUGMAPUBE1BU9BcDTSzyBNLwSoDmO6CZkLcFH9Bs84paOZml3laKVllABMETUFT0JxkC7F2HJqeWApWEcBsAErQFDSHOQXNHNDMIkcLCPdCYx98QTMlboq/oFnnFDSzyBNLweo/MKEJpPpY0BQ0hzkFzWWgmeWnkVSvsa5h9kWCpqCJoLkKaGaTJ5aCVQYwu4BUHwuaguYwp6C5k9DUCLMniiElaDb7Kmi2cwqa80Izy6jP0QLCfZKg2R1D0OzKKWjODc0lZWiE2R8NavgIml0xBM2unIJmKjSzyRNLweo/MCPAgaApaHbE6rIJmk1biNUFzVwMM/ekUrL6D0wiqAmagqaguRJoLi2foxSsAoAJuKA52a/ZJmh25RQ0p0Ezl3QNs08SNKf4NdsEza6cguaqoamfRvZFo0+GoDnRr9kmaHblFDS7oEmuVSo9sRSsMoAJguZYjC6/Zpug2ZVT0GxDM8t9GE+bjmtKvmKZM/pjV0FT0BQ0VwLNLPLEsoTM7L+Y2ffM7Ltm9j/qDRpD221mdtrMTpnZDZH9WjN7MrTdETZhJGzU+PvB/piZXTkrf++BCeFLLWgKmiBorgCaOWSsbYT5IHCNu/8K8H3gNgAzuxo4DLybahvwO8PGjQBfBI5S7WC7N7QDfBz4S3f/JeDzwOdmJS8CmCBoCpqC5tCnr9AceFJZRu7+P939XDj838BWqB8C7nH3N9z9BeA0sN/MdgGXuPuj7u7A3cBN0TnHQ/0PgOvq0ecklQHMIVQEzfF+T/NrtgmaXTkFzSzyOQpcZmYno3J0waz/ltFW3buBF6O2M8G2O9Tb9sY5AcJ/DbxzWsKUfcn7IXeqLW4dH1iF+mDDPWyhSgXNjeAX79tab49bH4cY5trCd1JcbeHbnXOUJM5Z7ha+uTTHnwy96u77JsYx+xPgFzqabnf3+4LP7cA54Kv1aR3+PsU+7ZyJKgOYERgFTUFT0CQzNDMpUyx3/8C0djM7AnwYuC5Ms6EaOe6J3LaAl4J9q8Men3PGzC4Afg54bVruMqbkQDwFBzQ9j+ONxejya7bN9ve5/JfKSYdtkl9nfE3PGyHmnJ7nknlaWSqH2UHgPwA3uvvfRU0ngMPhzvdVVDd3Hnf3l4HXzexAuD55M3BfdM6RUP914OEIwJ3qPzAngBEETUFT0Bz67DQ0fRh4dllOvw38LPCgmX3HzP4bgLs/DdwLPAM8ANzi7tvhnE8AX6K6EfQ8o+ueXwbeaWangU8Dt85KXsSUfNIUXNPzuE3Tc03PF5yek0fr+Nlj+BOgSW3HgGMd9pPANR32HwMfmSd//0eYQZNGk9PaNNKc5Nds00izK+f5NdJcVgZrmZLvtMoA5gwwTmsTNCf5NdsEza6c5wc0yQHN1Ol41rtM61cZwARBk+nwETTbNkGzEWIKNHNJI8y+KBGM09oEzUl+zTZBsyunoJkkTywFqwxgVrs1VXVBU9AUNLNCM5c0wuyDYsgJmoLmtLj1saC5fmg6sO1ppWD1H5gwelcFzYlxBU1Bs26bF5q5pBFmJDPbNLNvm9nXwvGlZvagmT0XHt8R+c61Lt10NSEoaE6OK2gKmnXbjkDTPa0UrHlGmJ8Eno2ObwUecve9wEPheNF16aarBUFBc3JcQVPQrNuSoZlJGmEGmdkW8CGqnxfVOsRoLbnjNNeYm3dduumK/2cSNAVNOmyC5hLQZHn5HKVgpY4wvwB8hhFyAK4IP2wnPF4e7IusS9eQmR2t18r7yXb4fb2g2eqHoJkUtz4WNKdCc1kZYNueVErWTGCa2YeBs+7+zcSYXWN8n2IfN7rf5e773H3fRZs/M1qlWdBs9UPQTIpbHwuandDMJXNPKiUrZYT5PuBGM/tz4B7g/Wb2FeCVMM0mPJ4N/ousSzdZQ4AImoJmhz8dNkFz/dD0OUrBmglMd7/N3bfc/UqqmzkPu/vHaK4ld4TmGnPzrks3qxOAoDnWJmgKmjmgmUU++l7OKgVrmb/D/CxwvZk9B1wfjllwXbrZEjQFzWkx6LAJmknQzCXztFKy5loP090fAR4J9R8C103wm2tduhlZq3c9XgNz4PiGjT4NUVu1/mXTZq71NLWeZkessfjn6XqauVT46DFF5fzSpz061EhTI80ufzpsk/w645+HI80ccnSXvFcSNMNzntImaAqaOwzNpFKwygDmFDAKmoKmoNlqmxOao7n+ctKfFfVBLVgJmqEuaAqaROe02+aGZgbV379ZpWD1H5jQhF99LGgKmikx6LAJmuPQXFZO9XlPKQWrHGAKmoLmojnpsAmaWQd7Rtp0XFPylasDgoKmoClo9g6aDAZppWAVAExGL7Kg2fIXNOfKSYdN0MwjR1PyXknQFDQFzV5DU1PyPsiJQCdoCpoZctJhO5+hmUv192xWKVj9BybgDUAKmoJmhpx02ATNJZQISwFzPRI0O2IJmoLmMOeS0FxWjnaN7I0ChATNjliCpqA5zLkgNDNJ1zB7oSYMBc2OWIKmoDnMuYPQ1JS8J/LwrRc0BU1BMz80c8ipPu8ppWCVAUwQNGN/QXOiv6BZ55wDmlnko++TRpg7KCcCkaApaM72FzTrnGuGpoDZEzUAJmgKmrP9Bc06ZwI0c8iB7UFaKVhlABMETUFT0EzIuXPQ9Op7mVIyyMz+vZm5mV0W2W4zs9NmdsrMbojs15rZk6HtjrAJI2Gjxt8P9sfM7MpZecsAZgxIQVPQnMNf0KxzToMmebSmKbmZ7aHaePEvItvVVLvavhs4CNxpZpuh+YvAUaodbPeGdoCPA3/p7r8EfB743KzcBQCzA5CCpqA5h7+gWedc5c8iWedd8s8Dn6H5DA4B97j7G+7+AtXOtPvNbBdwibs/6u4O3A3cFJ1zPNT/ALiuHn1OUv+BWb8RIGiG10PQFDRTcq4fmp5WlpCZ3Qj8X3d/otW0G3gxOj4TbLtDvW1vnOPu54C/Bt45Lf9c2+zulNy9+mJvWAVI26C5nW5ss8qf0DYYwMbG6I0absVrTZs72sK3o23YD7SFLwlx6+Ox+D3ewjeX0mF4mZmdjI7vcve76gMz+xPgFzrOux34j8C/7mjreiaTnmH0ak1s61QRwIT5oRlOEjQFTUFzFjRzyB22t1O9X3X3fZND+Qe67Gb2j4GrgCfCzHkL+JaZ7acaOe6J3LeAl4J9q8NOdM4ZM7sA+DngtWkd7/+UHIZTa3dPnp67u6bnXbGGtihWXdf0fHYMOmyT/Drj93B6nksrnpK7+5Pufrm7X+nuV1IB71fd/f8BJ4DD4c73VVQ3dx5395eB183sQLg+eTNwXwh5AjgS6r8OPByuc05UGcCEhaA59Bc0Bc2Jz32BnHTYBM2VA3N6an8auBd4BngAuMXd6yHvJ4AvUd0Ieh64P9i/DLzTzE4DnwZunZWnjCl5a2qta5qans+Kq+l5d85RkihWlguZnusOeHrGapQZHx8DjnX4nQSu6bD/GPjIPDn7P8IcG9lppKmRZlrcFH+NNDPJwX2QVEpW74Hp0ARcdCxoMvILL5agKWim5FwJNPXTyH6oAqOgKWg2+ypotnPuIDTr76i22e2HBE1a50bHgqagOcy5IDRzqP5OzCoFqwxgxqATNAXNqE3Q7Mo5JzQzyQeDpFKyCgBmE2CCJq1zo2NBU9Ac5pwDmlnko++CRpg7KKcbdIKmoBm1CZpdOdcITYc1Lr6xY+o/MAEfuKAZ9UPQ7I4haHblnA3NHHLAt7eTSskqApggaAqak/yabYJmV841QNPD93NNCwjvlMoAZniRBU0a/RA0u2MIml05p0Azk3zgSaVklQFMWBs0u2KcN9Bsn8ukL/hPITRn9OOnGpq5dB6MMG3G4hw7LjN7HTi10/1I1GXAqzvdiTlUUn9L6iuU1d+/7+7vWiaAmT1A9ZxT9Kq7H5zt1j+VAMyT09bO65NK6iuU1d+S+grl9VdKUzlTckmSpB2WgClJkpSoEoB512yX3qikvkJZ/S2pr1Bef6UE9f4apiRJUl9UwghTkiSpFxIwJUmSEtVbYJrZQTM7ZWanzWzm5kQr7MfvmNlZM3sqsl1qZg+a2XPh8R1R222hz6fM7IbIfq2ZPRna7gg72OXu6x4z+1Mze9bMnjazT/a1v2b2FjN73MyeCH39zb72NcqzaWbfNrOv9b2v0opUb0fbpwJsUu3u9g+Ai4AngKt3qC//AvhV4KnI9p+BW0P9VuBzoX516OvFVPsnPw9shrbHgX9G9UOM+4FfW0Ffd1FtOwrws8D3Q596198Q9+2hfiHwGHCgj32N+vxp4L8DX+vz50BldaWvI8z9wGl3/4G7/wS4Bzi0Ex1x9z9jfHP3Q8DxUD8O3BTZ73H3N9z9BaptPfeb2S7gEnd/1Ktvzd3ROTn7+rK7fyvUXweeBXb3sb9e6W/D4YWheB/7CmBmW8CHqLZrrdXLvkqrU1+BuRt4MTo+E2x90RVebRBPeLw82Cf1e3eot+0rk5ldCbyXauTWy/6GKe53gLPAg+7e274CXwA+w+gX/vS4r9KK1Fdgdl3XKeHvnyb1e63Px8zeDvwh8Cl3/5tprh22tfXX3bfd/T3AFtUIbGzv6Eg71lcz+zBw1t2/mXpKh23tnwMpv/oKzDPAnuh4C3hph/rSpVfC9IrweDbYJ/X7TKi37dllZhdSwfKr7v5Hfe8vgLv/FfAIcLCnfX0fcKOZ/TnV5aH3m9lXetpXaYXqKzC/Aew1s6vM7CLgMHBih/sU6wRwJNSPAPdF9sNmdrGZXQXsBR4P07XXzexAuCt6c3RONoXYXwaedfff6nN/zexdZvbzof5W4APA9/rYV3e/zd233P1Kqs/iw+7+sT72VVqxdvqu06QCfJDqLu/zwO072I/fA14G3qQaIXwceCfwEPBceLw08r899PkU0R1QYB/wVGj7bcKvrDL39Z9TTfG+C3wnlA/2sb/ArwDfDn19CvhPwd67vrb6/S8Z3SXvdV9V8hf9NFKSJClRfZ2SS5Ik9U4CpiRJUqIETEmSpEQJmJIkSYkSMCVJkhIlYEqSJCVKwJQkSUrU/wdEuPKLDN4MOAAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 2 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.imshow(uArr)\n",
    "plt.colorbar()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "5ba8359d",
   "metadata": {},
   "source": [
    "## 四、三维图"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "id": "f4a1ab9f",
   "metadata": {},
   "outputs": [],
   "source": [
    "A = np.identity(xArr.size - 2)\n",
    "B = np.roll(A, 1)\n",
    "C = np.roll(A, -1)\n",
    "B[0, 0] = 0\n",
    "C[-1, -1] = 0"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "id": "390a8f13",
   "metadata": {},
   "outputs": [],
   "source": [
    "A1 = b * A + c * B + a * C\n",
    "A2 = e * A + f * B + d * C\n",
    "A3 = np.linalg.inv(A1) @ A2"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "id": "9e331ed1",
   "metadata": {},
   "outputs": [],
   "source": [
    "uArr = np.zeros((xArr.size - 2, tArr.size))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "id": "5981d586",
   "metadata": {},
   "outputs": [],
   "source": [
    "uArr[:, 0] = uBegin[1:-1]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "id": "cbfee24c",
   "metadata": {},
   "outputs": [],
   "source": [
    "uCurrent = uArr[:, 0]\n",
    "for i in range(1, n):\n",
    "    uNext = A3 @ uCurrent\n",
    "    uArr[:, i] = uNext\n",
    "    uCurrent = uNext"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "id": "74d52321",
   "metadata": {},
   "outputs": [],
   "source": [
    "XX, YY = np.meshgrid(tArr, xArr[1:-1])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "id": "387760b7",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(2000, 5001)"
      ]
     },
     "execution_count": 23,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "XX.shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "id": "9aa0dc1b",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(2000, 5001)"
      ]
     },
     "execution_count": 24,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "uArr.shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "id": "7fc38eb8",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Using matplotlib backend: Qt5Agg\n"
     ]
    }
   ],
   "source": [
    "%matplotlib"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "id": "f8854171",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<mpl_toolkits.mplot3d.art3d.Poly3DCollection at 0x16197693100>"
      ]
     },
     "execution_count": 26,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "from mpl_toolkits.mplot3d import Axes3D\n",
    "fig = plt.figure(figsize = (12, 8))\n",
    "ax = fig.add_subplot(111, projection='3d')\n",
    "ax.plot_surface(XX, YY, uArr, cmap = plt.cm.hot)"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.9.7"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}
